/**
 * @file   EquationSolver.h
 * @author ZL <zuoling@Black>
 * @date   Thu Oct  6 14:28:47 2022
 * 
 * @brief  
 * 
 * 
 */

#include <iostream>  //std::cout
#include <limits>    //std::numeric_limits
#include <cmath>
#include <cfloat>

#define pi M_PI /// define const pi
const double eps = std::numeric_limits<double>::epsilon();/// define min number -> eps.
 //#define eps DBL_EPSILON

/** 
 * 
 * 
 * 
 * @return 创建函数的虚类
 */
class Function
{
public:
    virtual double operator()(double _x) = 0;
    virtual double diff(double _x){return 0;};/// 如果函数不需要导函数，默认导函数为0.
    
};

/** 
 * 
 * 
 * 
 * @return 创建求解器的虚类
 */
class EquationSolver
{
public:
    virtual double solve() = 0;
};

/** 
 * 
 * 
 * @param _x 
 * 
 * @return 返回输入数的正负
 */
int sgn(double _x)
{
    if (_x < 0)
	return -1;
    else if (_x = 0)
	return 0;
    else
	return 1;
    
};

/** 
 * 
 * 
 * @param _a 
 * @param _b 
 * @param _delta 
 * @param _M 
 * @param _f 
 * 
 * @return Bisection method 
 */
class Bisection : public EquationSolver
{
private:
    double a;
    double b;
    double delta;
    int M;
    Function &f;
public:
    Bisection(double _a,
	      double _b,
	      double _delta,
	      int _M,
	      Function &_f):
	a(_a),
	b(_b),
	delta(_delta),
	M(_M),
	f(_f){};
    
    double solve()
	{
	    double u = f(a);
	    double v = f(b);
	    double w;
	    double h;
	    double c;
	    for (int k = 0; k < M; k++)
	    {
		h = b-a;
		c = a+h/2;
		w = f(c);
		if (fabs(h) < delta || fabs(w) < eps)
		    break;
		else if (sgn(w) != sgn(u))
		{
		    b = c;
		    v = w;
		}
		else
		{
		    a = c;
		    u = w;
		}
	    }
	    return c;
	}
};

/** 
 * 
 * 
 * @param _x_0 
 * @param _delta 
 * @param _M 
 * @param _f 
 * 
 * @return Newton method
 */
class Newton : public EquationSolver
{
private:
    double x_0;
    double delta;
    int M;
    Function &f;
    
public:
    Newton(double _x_0, double _delta, int _M, Function &_f):
	x_0(_x_0),
	delta(_delta),
	M(_M),
	f(_f){};
    double solve()
	{
	    double x = x_0;
	    double u;
	    for (int k = 0; k <=M; k++)
	    {
		u = f(x);
		if (fabs(u) < eps)
		    break;
		x = x-u/f.diff(x);
	    }
	    return x;
	}
};

/** 
 * 
 * 
 * @param _x_0 
 * @param _x_1 
 * @param _delta 
 * @param _M 
 * @param _f 
 * 
 * @return Secant method 
 */
class Secant : public EquationSolver
{
private:
    double x_0;
    double x_1;
    double delta;
    int M;
    Function &f;
    
public:
    Secant(double _x_0,
	   double _x_1,
	   double _delta,
	   int _M,
	   Function &_f):
	x_0(_x_0),
	x_1(_x_1),
	delta(_delta),
	M(_M),
	f(_f){};
    double solve()
	{
	    double x_n = x_1;
	    double x_m = x_0;
	    double u = f(x_n);
	    double v = f(x_m);
	    double t;
	    for (int k = 2; k <= M; k++)
	    {
		if (fabs(u) > fabs(v))
		{
		    t = x_n; x_n = x_m; x_m = t;
		    t = u; u = v; v = t;
		}
		double s = (x_n-x_m)/(u-v);
		x_m = x_n;
		v = u;
		x_n = x_n-u*s;
		u = f(x_n);
		if (fabs(x_n-x_m) < delta || fabs(u) < eps)
		    break;
	    }
	    return x_n;
	}
};
